Things that should only have to be defined once, here at the top¶
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import igrins_mod as ig # Custom module file for igrins shenanigans
import glob
import os
import warnings
import astropy.units as u
from astroquery.nist import Nist # atomic lines
# from astroquery.linelists.cdms import CDMS # molecular lines?
from lmfit import Model, Parameters
from lmfit.models import VoigtModel, LorentzianModel, SkewedGaussianModel, SkewedVoigtModel
from muler.igrins import IGRINSSpectrum, IGRINSSpectrumList
from tqdm import tqdm
# from scipy.integrate import trapz, simpson
# from scipy.signal import savgol_filter
from astropy.io import fits
# Plotting Parameters
plt.rcParams['font.size'] = 20
plt.rcParams['axes.labelsize'] = 18
plt.rcParams['xtick.labelsize'] = 18
plt.rcParams['ytick.labelsize'] =18
plt.rcParams['legend.fontsize'] = 16
plt.rcParams['figure.titlesize'] = 20
plt.rcParams['axes.labelweight']='bold'
plt.rcParams['axes.linewidth'] = 3
plt.rcParams['xtick.major.size'] = 10
plt.rcParams['xtick.minor.visible'] = True
plt.rcParams['xtick.minor.size'] = 5
plt.rcParams['ytick.major.size'] = 10
plt.rcParams['ytick.minor.visible'] = True
plt.rcParams['ytick.minor.size'] = 5
plt.rcParams['ytick.direction'] = 'in'
plt.rcParams['xtick.direction'] = 'in'
%matplotlib inline
# Size of 1 spectral resolution element
# IGRINS Spectral Resolution
spec_res = 0.00001
# Reduced and order-merged data filepath
# Laptop Path
data_path = "C:\\Users\\Savio\\Documents\\GitHub\\IGRINS-Spectra\\IGRINS_Merged"
# File path for figures to live in
fig_path = "C:\\Users\\Savio\\Documents\\GitHub\\IGRINS-SpectraIGRINS_figs\\standards_spectra"
# Create the folder if it doesn't exist
if not os.path.exists(fig_path):
os.makedirs(fig_path)
# Nicole's merged K-band spectra of some Taurus Standards
# merged_standard_files = glob.glob(data_path + "/merged_standards/m*.fits")
standard_table = pd.read_csv('C:\\Users\\Savio\\Documents\\IGRINS-Spectra\\standard_table_v3.txt', index_col=0) # csv of standards with file and Spectral Type, c/v TBA
proto_table = pd.read_csv('C:\\Users\\Savio\\Documents\\IGRINS-Spectra\\protostar_table.txt', index_col=0)
# mask for spectra with well defined peaks and little-to-no line blending
feature_mask = standard_table['Name'].isin( ['LkCa19', 'HBC427', 'LkCa4','Hubble4','V830Tau', 'L1551-55','L1551-51', 'Anon1','JH108', 'UXTauB', 'LkCa1'] )
standard_table = standard_table[feature_mask].reset_index(drop=True)
standards_path = standard_table['File']
standard_list = standard_table['File'].values
proto_path = proto_table['File']
proto_list = proto_table['File'].values
standard_stack = [] # initialize list of IGRINSSpectrum objects
proto_stack = []
num_orders = 25
order_len = 2048
for i in tqdm(range(len(standard_table))):
spec = IGRINSSpectrumList.read(standards_path[i])\
.trim_overlap(pivot=0.5)\
.stitch()
standard_stack.append(spec)
for i in tqdm(range(len(proto_table))):
spec = IGRINSSpectrumList.read(proto_path[i])\
.trim_overlap(pivot=0.5)\
.stitch()
proto_stack.append(spec)
# Determine the maximum length of flux arrays for the standards
max_flux_length = max(len(spec.flux.value) for spec in standard_stack)
max_wavelen_length = max(len(spec.wavelength.value) for spec in standard_stack)
max_snr_length = max(len(spec.snr.value) for spec in standard_stack)
max_unc_length = max(len(spec.uncertainty.array) for spec in standard_stack)
# Initialize stacks with NaN values
wavelen_stack = np.full((max_wavelen_length, len(standard_table)), np.nan)
flux_stack = np.full((max_flux_length, len(standard_table)), np.nan)
snr_stack = np.full((max_snr_length, len(standard_table)), np.nan)
flux_unc_stack = np.full((max_unc_length, len(standard_table)), np.nan)
for i in tqdm(range(len(standard_table))):
# convert angstrom to micron
wavelen = standard_stack[i].wavelength.value*1e-4
flux = standard_stack[i].flux.value
snr = standard_stack[i].snr.value
raw_flux_err = standard_stack[i].uncertainty.array
wavelen_mask = (wavelen > 2.1) & (wavelen < 2.39)
flux_mask = (flux > 0)
# flux_err_mask = (raw_flux_err < 0.1)
mask = (wavelen_mask) & (flux_mask) #& (flux_err_mask)
wavelen = wavelen[mask]
flux = flux[mask]
snr = snr[mask]
raw_flux_err = raw_flux_err[mask]
wavelen_stack[:len(wavelen),i] = wavelen
flux_stack[:len(flux),i] = flux
snr_stack[:len(snr),i] = snr
flux_unc_stack[:len(raw_flux_err),i] = raw_flux_err
# Directly query NIST to find line features in K-band
with warnings.catch_warnings(): # Ignore warnings
warnings.simplefilter('ignore')
lines_table = Nist.query(2.08*u.um,2.35*u.um,
linename = 'Na I, Sc I, Si I, Fe I, Al I, Mg I, Ca I, H I, Ti I',
energy_level_unit='eV',output_order='wavelength')
igrins_wav_cut = (lines_table['Observed'] > 2.08) & (lines_table['Observed'] < 2.35)
lines_table = lines_table[igrins_wav_cut]
# lines_table = pd.read_csv('lines_table.txt')
# Make masks for the table of all the lines just in case I want to peek at certain transitions/wavelengths
na1_mask = lines_table['Spectrum'] == 'Na I'
sc1_mask = lines_table['Spectrum'] == 'Sc I'
si1_mask = lines_table['Spectrum'] == 'Si I'
fe1_mask = lines_table['Spectrum'] == 'Fe I'
al1_mask = lines_table['Spectrum'] == 'Al I'
mg1_mask = lines_table['Spectrum'] == 'Mg I'
ca1_mask = lines_table['Spectrum'] == 'Ca I'
h1_mask = lines_table['Spectrum'] == 'H I'
ti1_mask = lines_table['Spectrum'] == 'Ti I'
# Just add all the masks to a list for the sake of my plotting a few cells down
mask_list = [na1_mask,sc1_mask,si1_mask,fe1_mask,al1_mask,mg1_mask,ca1_mask,h1_mask,ti1_mask]
color_list = ['purple', 'orange', 'green', 'blue', 'brown', 'crimson', 'olive', 'cyan', 'darkgreen']
na_reg_lines = sorted([lines_table[na1_mask][0]['Observed'],
lines_table[sc1_mask][19]['Observed'],
lines_table[si1_mask][2]['Observed'],
lines_table[sc1_mask][20]['Observed'],
lines_table[na1_mask][1]['Observed']])
na_reg_lines = np.array(na_reg_lines)
ca_reg_lines = list(lines_table[ca1_mask]['Observed'][:]) # list of Ca I lab wavelengths
ca_reg_lines.append(lines_table[fe1_mask]['Observed'][104])
ca_reg_lines = np.array(sorted(ca_reg_lines))
100%|██████████| 11/11 [00:02<00:00, 4.22it/s] 100%|██████████| 48/48 [00:10<00:00, 4.53it/s] 100%|██████████| 11/11 [00:00<00:00, 843.83it/s]
ca_reg_lines = list(lines_table[ca1_mask]['Observed'][:]) # list of Ca I lab wavelengths
ca_reg_lines.append(lines_table[fe1_mask]['Observed'][104])
ca_reg_lines = np.array(sorted(ca_reg_lines))
# Define the region for fitting
line_name = lines_table[ca1_mask][0]['Spectrum'] # Species
line_center = ca_reg_lines[1] # Wavelength
# from igrins_mod import local_continuum_fit
continuum_stack = []
continuum_err_stack = []
norm_flux_stack = []
norm_flux_err_stack = []
reg_idx_stack = []
# regions [left point,width]
regions = [(-1000,11), (-500,250), (190,51), (500,11) , (1000,11)]
# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
poly_deg = 3
for i in tqdm(range(len(standard_table))):
fig = plt.figure(figsize=(15,3))
# Perform the continuum fit
continuum, continuum_err, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
flux_stack[:,i],
flux_unc_stack[:,i],
poly_order = poly_deg,
line_center = line_center,
spec_res = spec_res,
regions = regions)
# Append indices to list of indices
reg_idx_stack.append(region_indices)
# Append to list of local continuum arrays
continuum_stack.append(continuum)
continuum_err_stack.append(continuum_err)
# Normalize flux by dividing the original flux by the continuum
norm_flux = flux_stack[region_indices[0][0]:region_indices[n][1], i] / continuum
norm_flux_stack.append(norm_flux)
# Plot the raw flux and continuum
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1], i], flux_stack[region_indices[0][0]:region_indices[n][1], i])
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1], i], continuum_stack[i])
# Highlight the regions and Ca I lines
for j in range(len(region_indices)):
plt.axvspan(wavelen_stack[region_indices[n - j][0], i], wavelen_stack[region_indices[n - j][1], i], color='red', alpha=0.2)
for k in range(len(ca_reg_lines)):
plt.axvline(ca_reg_lines[k], ymin=0.8, ymax=0.9, color='black')
plt.title(f"{line_name} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
plt.show()
0%| | 0/11 [00:00<?, ?it/s]
9%|▉ | 1/11 [00:00<00:03, 2.52it/s]
18%|█▊ | 2/11 [00:00<00:04, 2.17it/s]
27%|██▋ | 3/11 [00:01<00:02, 2.79it/s]
36%|███▋ | 4/11 [00:01<00:02, 3.15it/s]
45%|████▌ | 5/11 [00:01<00:01, 3.26it/s]
55%|█████▍ | 6/11 [00:01<00:01, 3.42it/s]
64%|██████▎ | 7/11 [00:02<00:01, 3.65it/s]
73%|███████▎ | 8/11 [00:02<00:00, 3.60it/s]
82%|████████▏ | 9/11 [00:02<00:00, 3.47it/s]
91%|█████████ | 10/11 [00:03<00:00, 3.64it/s]
100%|██████████| 11/11 [00:03<00:00, 3.01it/s] 100%|██████████| 11/11 [00:03<00:00, 3.01it/s]
7/31/24¶
compare methods of fitting (letting things float vs fixing them)
compare fitted delta_lam vs lab delta_lam
consider printing offsets in centers from my original fitting method where centers are free and the new method where we fix the centers to be equal to lab wavelengths
lab wavelengths might be slightly off hence the need for all this
fit without beta, calculate the deltas using a high S/N line as the reference then fiddle with the lab wavelengths afterwards
8/5¶
test the delta lam for the not blended lines
Compare offset values between lab and fit values for the spectra with good composite fits so pick the *good spectra* and average the deltas of the deltas
look at the fractional sigmas of the fits (sigma/sigma_avg) where sigma is free parameter
should we really assume sigmas should be the same for all lines???
Method 3¶
ca_reg_lines
array([2.26141 , 2.262601, 2.26311 , 2.265741])
# Initialize storage lists
params_stack = []
beta_stack = []
params_error_stack = []
best_model_stack = []
result_stack = []
dely_stack = []
# For testing purposes
mean_sigma_stack = []
# ca_obs = ca_lines*(1+standard_table['v/c'][0]) # observed wavelengths (theoretically)
for i in range(len(standard_table)):
# Define initial parameters for Gaussian fitting
params = Parameters()
# Define Model
# model = Model(ig.four_gaussian_beta, nan_policy='omit')
gauss_model = Model(ig.four_gaussian_beta, nan_policy='omit')
# extra_model = SkewedVoigtModel(prefix='m4_')
model = gauss_model
params.add('amp1', value = -2e-5, max = 0)
params.add('c1', value = ca_reg_lines[0], max=ca_reg_lines[0]+10*spec_res)
params.add('std1', min=0, value = 1.5e-4)
params.add('amp2', value = -2e-5, max = 0)
params.add('c2', value = ca_reg_lines[1], max=ca_reg_lines[1]+10*spec_res)
params.add('std2', min=0, value = 1.5e-4)
params.add('amp3', value = -2e-5, max = 0)
params.add('c3', value = ca_reg_lines[2], max=ca_reg_lines[2]+10*spec_res)
params.add('std3', min=0, value = 1.5e-4)
params.add('amp4', value = -2e-5, max = 0)
params.add('c4', value = ca_reg_lines[3], max=ca_reg_lines[3]+10*spec_res)
params.add('std4', value = 1.5e-4)
# hold widths to be the same
# params['std2'].expr = 'std1'
# params['std3'].expr = 'std1'
# params['std4'].expr = 'std1'
# Add parameters for the Voigt profile
# params.add('m4_amplitude', value=-5e-5, max=0)
# params.add('m4_center', value=ca_reg_lines[3])
# params.add('m4_sigma', value=1.5e-4)
# params.add('m4_gamma', value=1.5e-4)
# params.add('m4_skew', value =1.5e-4)
params.add('beta', value=1.) # beta = (1+v/c)
# print(f"Length of norm_flux_stack[{i}]: {len(norm_flux_stack[i])}")
# print(f"Length of wavelen_stack slice: {len(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i])}")
result = model.fit(norm_flux_stack[i], params,
x = wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i])
# weights = 1/flux_unc_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i])
result_stack.append(result) # lmfit ModelResult objects
dely = result.eval_uncertainty(sigma=3) # lmfit uncertainty
dely_stack.append(dely)
# Append the best_model, error and parameters stacks
# could probably make a loop to the tune of for value in result etc...
params_stack.extend([result.params['amp1'].value, result.params['c1'].value, result.params['std1'].value,
result.params['amp2'].value, result.params['c2'].value, result.params['std2'].value,
result.params['amp3'].value, result.params['c3'].value, result.params['std3'].value,
result.params['amp4'].value, result.params['c4'].value, result.params['std4'].value])
# result.params['m4_amplitude'].value, result.params['m4_center'].value, result.params['m4_sigma'].value])
beta_stack.append(result.params['beta'].value)
params_error_stack.extend([result.params['amp1'].stderr, result.params['c1'].stderr, result.params['std1'].stderr,
result.params['amp2'].stderr, result.params['c2'].stderr, result.params['std2'].stderr,
result.params['amp3'].stderr, result.params['c3'].stderr, result.params['std3'].stderr,
result.params['amp4'].stderr, result.params['c4'].stderr, result.params['std4'].stderr])
# result.params['m4_amplitude'].stderr, result.params['m4_center'].stderr, result.params['m4_sigma'].stderr])
best_model_stack.append(result.best_fit)
mean_sigma_stack.append(np.mean(a=(result.params['std1'].value, result.params['std2'].value, result.params['std3'].value, result.params['std4'].value)))
# print(result.fit_report())
# print(f"Widths: {result.params['std1'].value:.8f}, {result.params['std2'].value:.8f}, {result.params['std3'].value:.8f}, {result.params['std4'].value:.8f}")
# print(f"Mean Width: {np.mean(a= (result.params['std1'].value, result.params['std2'].value, result.params['std3'].value, result.params['std4'].value)):.8f}")
params_arr = np.array(params_stack)
# 0: amplitude
# 1: centers
# 2: sigmas
amps = params_arr[0::3] # grab every amplitude
amps_err = params_error_stack[0::3]
centers = params_arr[1::3] # grab every center wavelength
center_err = params_error_stack[1::3]
sigmas = params_arr[2::3] # grab every sigma
sigmas_err = params_error_stack[2::3]
# print(centers)
# Gaussian model fits for each source
ca1_2_2614_fits = []
ca1_2_2631_fits = []
ca1_2_2657_fits = []
fe1_2_2626_fits = []
num_gauss = 4
# # Getting the individual Gaussians into their own lists
for i in range(len(standard_table)):
ca1_2_2614_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[0::num_gauss][i], centers[0::num_gauss][i]*beta_stack[i], sigmas[0::num_gauss][i]))
fe1_2_2626_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[1::num_gauss][i], centers[1::num_gauss][i]*beta_stack[i], sigmas[1::num_gauss][i]))
ca1_2_2631_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[2::num_gauss][i], centers[2::num_gauss][i]*beta_stack[i], sigmas[2::num_gauss][i]))
ca1_2_2657_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[3::num_gauss][i], centers[3::num_gauss][i]*beta_stack[i], sigmas[3::num_gauss][i]))
# Plotting each component Gaussian and the best model fit
for i in range(len(standard_table)):
# Create subplots with adjusted size ratios
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], best_model_stack[i], c='red', label=r'Best Fit $\pm 3 \sigma$')
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], norm_flux_stack[i], c='black', alpha=0.5)
# ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], result_stack[i].init_fit, '--', c='blue' , alpha=0.2, label='initial fit')
ax1.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
best_model_stack[i]-dely_stack[i], best_model_stack[i]+dely_stack[i],
alpha=0.2, color = 'red' )
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2614_fits[i], ls='-.',
label=r'$\bf{Ca~I~2.2614~\mu m}$' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[0::num_gauss][i]*beta_stack[i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")
# rf"$\pm {amps_err[0::num_gauss][i]:.3e}$"
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], fe1_2_2626_fits[i], ls='-.',
label=r'$\bf{Fe~I~2.2626~\mu m}$' "\n" rf"A={amps[1::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[1::num_gauss][i]*beta_stack[i]:.5f}" "\n" rf"$\sigma$= {sigmas[1::num_gauss][i]:.3e}")
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2631_fits[i], ls='-.',
label=r'$\bf{Ca~I~2.2631~\mu m}$' "\n" rf"A={amps[2::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[2::num_gauss][i]*beta_stack[i]:.5f}" "\n" rf"$\sigma$= {sigmas[2::num_gauss][i]:.3e}")
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2657_fits[i], ls='-.', color='goldenrod',
label=r'$\bf{Ca~I~2.2657~\mu m}$' "\n" rf"A={amps[3::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[3::num_gauss][i]*beta_stack[i]:.5f}" "\n" rf"$\sigma$= {sigmas[3::num_gauss][i]:.3e}")
for j in range(len(ca_reg_lines)):
ax1.axvline(ca_reg_lines[j], ymin=0.05, ymax=0.95,ls='--', color='black')
plt.suptitle(rf"Ca Region {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
ax1.set_ylabel('Flux + Constant')
ax1.legend(bbox_to_anchor=(1,1))
ax1.set_title(r'Free $\sigma$ and yes $\beta$ offset')
# Calculate residuals
residuals = (norm_flux_stack[i] - best_model_stack[i])
# sqsum_res = np.sum(residuals**2)
# sqsum = np.sum((norm_flux_stack[i]-np.mean(norm_flux_stack[i])**2))
# R2 = 1-(sqsum_res/sqsum)
# Plot residuals
# ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
result_stack[i].plot_residuals(ax=ax2, datafmt='r.')
ax2.axhline(y=0, color='grey', linestyle='--')
# ax2.set_ylim(-0.1,0.1)
ax2.set_xlabel(r'Wavelength [$\mu$m]')
ax2.set_ylabel('res')
ax2.set_title(rf"$R^2$={result_stack[i].rsquared:.3f}, " r'$\beta \equiv 1+\frac{v}{c}= $' f"{beta_stack[i]:.8f}", fontsize=12)
plt.show()
delt_diff_stack = []
ref_line = 0 # reference spectral feature index
for i in range(len(standard_table)):
# To store the differences for each spectrum
delta_diff_spectrum = []
for j in range(len(ca_reg_lines)):
# fit line center and ref line center offset
fit_ref_diff = np.abs(centers[j::num_gauss][i] - centers[ref_line::num_gauss][i])
# print(fit_ref_diff)
# lab line center and ref line center offset
lab_ref_diff = np.abs(ca_reg_lines - ca_reg_lines[ref_line])[j]
# print(lab_ref_diff)
# difference between the delta wavelengths
delta_diff = np.abs(fit_ref_diff - lab_ref_diff)
# print(delta_diff)
# append to list of differences for spectrum i
delta_diff_spectrum.append(delta_diff)
# append the differences for this spectrum to the main list
delt_diff_stack.append(delta_diff_spectrum)
# creaty empty array with shape (num_of_sources, num_of_lines)
line_center_arr = np.empty(shape=(len(standard_table), len(ca_reg_lines)))
for i in range(len(standard_table)):
# take lab wavelengths from NIST then add the delta of the wavelength offsets
line_center_arr[i] = (ca_reg_lines + delt_diff_stack[i])
# create empty array with length of # of lines in spectral region
ca_reg_tweak = np.empty(shape=(len(ca_reg_lines)))
for j in range(len(ca_reg_lines)):
# average the wavelength for each line
# avg of line 1, line 2, etc.
ca_reg_tweak[j] = np.mean(line_center_arr[:,j])
print("Tweaked Wavelength Array:", line_center_arr)
print("avg:", ca_reg_tweak)
Tweaked Wavelength Array: [[2.26141 2.26260792 2.26315397 2.26575502] [2.26141 2.26263035 2.26314837 2.26583554] [2.26141 2.26262104 2.26318359 2.26579992] [2.26141 2.26264115 2.26317282 2.26579946] [2.26141 2.26265675 2.26316419 2.26581989] [2.26141 2.26262004 2.26315178 2.26579528] [2.26141 2.26260694 2.26316863 2.26580015] [2.26141 2.2626218 2.26316351 2.2658521 ] [2.26141 2.26262779 2.26316132 2.26580026] [2.26141 2.26262783 2.26315964 2.26579952] [2.26141 2.26265746 2.26318768 2.26577629]] avg: [2.26141 2.26262901 2.26316504 2.26580304]
ca_reg_tweak
array([2.26141 , 2.26262901, 2.26316504, 2.26580304])
Tweak Wavlengths¶
# Initialize storage lists
params_stack = []
beta_stack = []
params_error_stack = []
best_model_stack = []
result_stack = []
dely_stack = []
# For testing purposes
mean_sigma_stack = []
# ca_obs = ca_lines*(1+standard_table['v/c'][0]) # observed wavelengths (theoretically)
for i in range(len(standard_table)):
# adjust lab wavelengths by the average
# ca_reg_tweak = ca_reg_lines + np.mean(delt_diff_stack[i])
# Define initial parameters for Gaussian fitting
params = Parameters()
# Define Model
model = Model(ig.four_gaussian_beta, nan_policy='omit')
params.add('amp1', value = -2e-5, max = 0)
params.add('c1', value = ca_reg_tweak[0])
params.add('std1', min=0, value = 1.5e-4)
params.add('amp2', value = -2e-5, max = 0)
params.add('c2', value = ca_reg_tweak[1])
params.add('std2', min=0, value = 1.5e-4)
params.add('amp3', value = -2e-5, max = 0)
params.add('c3', value = ca_reg_tweak[2])
params.add('std3', min=0, value = 1.5e-4)
params.add('amp4', value = -2e-5, max = 0)
params.add('c4', value = ca_reg_tweak[3])
params.add('std4', value = 1.5e-4)
# tweak lab wavelengths
params['c1'].expr = f'{ca_reg_tweak[0]}'
params['c2'].expr = f'{ca_reg_tweak[1]}'
params['c3'].expr = f'{ca_reg_tweak[2]}'
params['c4'].expr = f'{ca_reg_tweak[3]}'
params.add('beta', value=1.) # beta = (1+v/c)
params['std2'].expr = 'std1' # constrain widths to be the same
params['std3'].expr = 'std1'
params['std4'].expr = 'std1'
result = model.fit(norm_flux_stack[i], params,
x = wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i])
# weights = 1/flux_unc_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i])
result_stack.append(result) # lmfit ModelResult objects
dely = result.eval_uncertainty(sigma=3) # lmfit uncertainty
dely_stack.append(dely)
# Append the best_model, error and parameters stacks
# could probably make a loop to the tune of for value in result etc...
params_stack.extend([result.params['amp1'].value, result.params['c1'].value, result.params['std1'].value,
result.params['amp2'].value, result.params['c2'].value, result.params['std2'].value,
result.params['amp3'].value, result.params['c3'].value, result.params['std3'].value,
result.params['amp4'].value, result.params['c4'].value, result.params['std4'].value])
beta_stack.append(result.params['beta'].value)
params_error_stack.extend([result.params['amp1'].stderr, result.params['c1'].stderr, result.params['std1'].stderr,
result.params['amp2'].stderr, result.params['c2'].stderr, result.params['std2'].stderr,
result.params['amp3'].stderr, result.params['c3'].stderr, result.params['std3'].stderr,
result.params['amp4'].stderr, result.params['c4'].stderr, result.params['std4'].stderr])
best_model_stack.append(result.best_fit)
mean_sigma_stack.append(np.mean(a=(result.params['std1'].value, result.params['std2'].value, result.params['std3'].value, result.params['std4'].value)))
# print(result.fit_report())
# print(f"Widths: {result.params['std1'].value:.8f}, {result.params['std2'].value:.8f}, {result.params['std3'].value:.8f}, {result.params['std4'].value:.8f}")
# print(f"Mean Width: {np.mean(a= (result.params['std1'].value, result.params['std2'].value, result.params['std3'].value, result.params['std4'].value)):.8f}")
params_arr = np.array(params_stack)
# 0: amplitude
# 1: centers
# 2: sigmas
amps = params_arr[0::3] # grab every amplitude
amps_err = params_error_stack[0::3]
centers = params_arr[1::3] # grab every center wavelength
center_err = params_error_stack[1::3]
sigmas = params_arr[2::3] # grab every sigma
sigmas_err = params_error_stack[2::3]
# print(centers)
# Gaussian model fits for each source
ca1_2_2614_fits = []
ca1_2_2631_fits = []
# ca1_2_2657_fits = []
fe1_2_2626_fits = []
num_gauss = 4
# # Getting the individual Gaussians into their own lists
for i in range(len(standard_table)):
ca1_2_2614_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[0::num_gauss][i], centers[0::num_gauss][i]*beta_stack[i], sigmas[0::num_gauss][i]))
fe1_2_2626_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[1::num_gauss][i], centers[1::num_gauss][i]*beta_stack[i], sigmas[1::num_gauss][i]))
ca1_2_2631_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[2::num_gauss][i], centers[2::num_gauss][i]*beta_stack[i], sigmas[2::num_gauss][i]))
ca1_2_2657_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[3::num_gauss][i], centers[3::num_gauss][i]*beta_stack[i], sigmas[3::num_gauss][i]))
# Plotting each component Gaussian and the best model fit
for i in range(len(standard_table)):
# ca_reg_tweak = ca_reg_lines + np.mean(delt_diff_stack[i])
# Create subplots with adjusted size ratios
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], best_model_stack[i], c='red', label=r'Best Fit $\pm 3 \sigma$')
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], norm_flux_stack[i], c='black', alpha=0.5)
# ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], result_stack[i].init_fit, '--', c='blue' , alpha=0.2, label='initial fit')
ax1.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
best_model_stack[i]-dely_stack[i], best_model_stack[i]+dely_stack[i],
alpha=0.2, color = 'red' )
# plot each componenet
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2614_fits[i], ls='-.',
label=r'$\bf{Ca~I~2.2614~\mu m}$' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu \cdot \beta$={centers[0::num_gauss][i]*beta_stack[i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")
# rf"$\pm {amps_err[0::num_gauss][i]:.3e}$"
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], fe1_2_2626_fits[i], ls='-.',
label=r'$\bf{Fe~I~2.2626~\mu m}$' "\n" rf"A={amps[1::num_gauss][i]:.3e}" "\n" rf"$\mu \cdot \beta$={centers[1::num_gauss][i]*beta_stack[i]:.5f}" "\n" rf"$\sigma$= {sigmas[1::num_gauss][i]:.3e}")
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2631_fits[i], ls='-.',
label=r'$\bf{Ca~I~2.2631~\mu m}$' "\n" rf"A={amps[2::num_gauss][i]:.3e}" "\n" rf"$\mu \cdot \beta$={centers[2::num_gauss][i]*beta_stack[i]:.5f}" "\n" rf"$\sigma$= {sigmas[2::num_gauss][i]:.3e}")
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2657_fits[i], ls='-.', color='goldenrod',
label=r'$\bf{Ca~I~2.2657~\mu m}$' "\n" rf"A={amps[3::num_gauss][i]:.3e}" "\n" rf"$\mu \cdot \beta$={centers[3::num_gauss][i]*beta_stack[i]:.5f}" "\n" rf"$\sigma$= {sigmas[3::num_gauss][i]:.3e}")
for j in range(len(ca_reg_lines)):
ax1.axvline(ca_reg_lines[j], ymin=0.05, ymax=1., color='black', ls='--', alpha=0.5)
ax1.axvline(ca_reg_tweak[j]*beta_stack[i], ymin=0.05, ymax=1., color='black', ls='--')
print(f"Offset from tweaked lab wavelength to fit center for line {j+1}: {np.abs(ca_reg_lines[j]*beta_stack[i]-(ca_reg_tweak[j])*beta_stack[i]):.6e}")
plt.suptitle(rf"Ca Region {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
ax1.set_ylabel('Flux + Constant')
ax1.legend(bbox_to_anchor=(1,1))
ax1.set_title(r'Free $\sigma$ and yes $\beta$ offset')
# Calculate residuals
residuals = (norm_flux_stack[i] - best_model_stack[i])
# sqsum_res = np.sum(residuals**2)
# sqsum = np.sum((norm_flux_stack[i]-np.mean(norm_flux_stack[i])**2))
# R2 = 1-(sqsum_res/sqsum)
# Plot residuals
# ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
result_stack[i].plot_residuals(ax=ax2, datafmt='r.')
ax2.axhline(y=0, color='grey', linestyle='--')
# ax2.set_ylim(-0.1,0.1)
ax2.set_xlabel(r'Wavelength [$\mu$m]')
ax2.set_ylabel('res')
ax2.set_title(rf"$R^2$={result_stack[i].rsquared:.3f}, " r'$\beta \equiv 1+\frac{v}{c}= $' f"{beta_stack[i]:.8f}", fontsize=12)
plt.show()
Offset from tweaked lab wavelength to fit center for line 1: 4.440892e-16 Offset from tweaked lab wavelength to fit center for line 2: 2.800799e-05 Offset from tweaked lab wavelength to fit center for line 3: 5.504385e-05 Offset from tweaked lab wavelength to fit center for line 4: 6.203871e-05
Offset from tweaked lab wavelength to fit center for line 1: 4.440892e-16 Offset from tweaked lab wavelength to fit center for line 2: 2.800787e-05 Offset from tweaked lab wavelength to fit center for line 3: 5.504361e-05 Offset from tweaked lab wavelength to fit center for line 4: 6.203844e-05
Offset from tweaked lab wavelength to fit center for line 1: 4.440892e-16 Offset from tweaked lab wavelength to fit center for line 2: 2.800880e-05 Offset from tweaked lab wavelength to fit center for line 3: 5.504544e-05 Offset from tweaked lab wavelength to fit center for line 4: 6.204050e-05
Offset from tweaked lab wavelength to fit center for line 1: 4.440892e-16 Offset from tweaked lab wavelength to fit center for line 2: 2.800757e-05 Offset from tweaked lab wavelength to fit center for line 3: 5.504303e-05 Offset from tweaked lab wavelength to fit center for line 4: 6.203778e-05
Offset from tweaked lab wavelength to fit center for line 1: 4.440892e-16 Offset from tweaked lab wavelength to fit center for line 2: 2.800700e-05 Offset from tweaked lab wavelength to fit center for line 3: 5.504191e-05 Offset from tweaked lab wavelength to fit center for line 4: 6.203652e-05
Offset from tweaked lab wavelength to fit center for line 1: 4.440892e-16 Offset from tweaked lab wavelength to fit center for line 2: 2.800815e-05 Offset from tweaked lab wavelength to fit center for line 3: 5.504417e-05 Offset from tweaked lab wavelength to fit center for line 4: 6.203907e-05
Offset from tweaked lab wavelength to fit center for line 1: 4.440892e-16 Offset from tweaked lab wavelength to fit center for line 2: 2.800966e-05 Offset from tweaked lab wavelength to fit center for line 3: 5.504712e-05 Offset from tweaked lab wavelength to fit center for line 4: 6.204240e-05
Offset from tweaked lab wavelength to fit center for line 1: 4.440892e-16 Offset from tweaked lab wavelength to fit center for line 2: 2.800915e-05 Offset from tweaked lab wavelength to fit center for line 3: 5.504612e-05 Offset from tweaked lab wavelength to fit center for line 4: 6.204127e-05
Offset from tweaked lab wavelength to fit center for line 1: 4.440892e-16 Offset from tweaked lab wavelength to fit center for line 2: 2.801078e-05 Offset from tweaked lab wavelength to fit center for line 3: 5.504934e-05 Offset from tweaked lab wavelength to fit center for line 4: 6.204490e-05
Offset from tweaked lab wavelength to fit center for line 1: 4.440892e-16 Offset from tweaked lab wavelength to fit center for line 2: 2.800722e-05 Offset from tweaked lab wavelength to fit center for line 3: 5.504234e-05 Offset from tweaked lab wavelength to fit center for line 4: 6.203701e-05
Offset from tweaked lab wavelength to fit center for line 1: 4.440892e-16 Offset from tweaked lab wavelength to fit center for line 2: 2.800887e-05 Offset from tweaked lab wavelength to fit center for line 3: 5.504559e-05 Offset from tweaked lab wavelength to fit center for line 4: 6.204067e-05
ew1_stack = []
ew2_stack = []
ew3_stack = []
ew4_stack = []
for i in range(len(standard_table)):
ew1 = np.trapz(np.abs(1-ca1_2_2614_fits[i]), wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew2 = np.trapz(np.abs(1-fe1_2_2626_fits[i]), wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew3 = np.trapz(np.abs(1-ca1_2_2631_fits[i]), wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
# ew4 = np.trapz(np.abs(1-ca1_2_2657_fits[i]), wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew1_stack.append(ew1)
ew2_stack.append(ew2)
ew3_stack.append(ew3)
# ew4_stack.append(ew4)
fig = plt.figure(figsize=(15,5))
# for i in range(len(standard_table)):
plt.plot(standard_table['Spectral_Type'], ew1_stack, ls='--', marker='o', c ='k', label = 'Ca I 2.2614')
plt.plot(standard_table['Spectral_Type'], ew2_stack, ls='--', marker='o', c ='red', label = 'Fe I 2.2626')
plt.plot(standard_table['Spectral_Type'], ew3_stack, ls='--', marker='o', c ='blue', label = 'Ca I 2.2631')
# plt.plot(standard_table['Spectral_Type'], ew4_stack, ls='--', marker='o', c ='green', label = 'Ca I 2.2657')
plt.xlabel('Source Spectral Type')
plt.ylabel(r'Equivalent Width [$\mu m$]')
plt.legend(bbox_to_anchor=(1,1))
plt.show()
fig = plt.figure(figsize=(15,5))
# for i in range(len(standard_table)):
plt.plot(standard_table['Spectral_Type'], ew1_stack, ls='--', marker='o', c ='k', label = 'Ca I 2.2614')
plt.plot(standard_table['Spectral_Type'], ew2_stack, ls='--', marker='o', c ='red', label = 'Fe I 2.2626')
plt.plot(standard_table['Spectral_Type'], ew3_stack, ls='--', marker='o', c ='blue', label = 'Ca I 2.2631')
# plt.plot(standard_table['Spectral_Type'], ew4_stack, ls='--', marker='o', c ='green', label = 'Ca I 2.2657')
plt.xlabel('Source Spectral Type')
plt.ylabel(r'Equivalent Width Ratios')
plt.legend(bbox_to_anchor=(1,1))
plt.show()
Try different models on 3rd Ca I line¶
# # Initialize storage lists
# params_stack = []
# beta_stack = []
# params_error_stack = []
# best_model_stack = []
# result_stack = []
# dely_stack = []
# # For testing purposes
# mean_sigma_stack = []
# # ca_obs = ca_lines*(1+standard_table['v/c'][0]) # observed wavelengths (theoretically)
# for i in range(len(standard_table)):
# # Define Model
# gauss_model = Model(ig.three_gaussian_beta, nan_policy='omit')
# # extra_model = SkewedVoigtModel(prefix='m4_')
# model = gauss_model
# # Define initial parameters for Gaussian fitting
# params = Parameters()
# params.add('amp1', value = -3e-5, max = 0)
# params.add('c1', value = ca_reg_tweak[0])
# params.add('std1', min=0, value = 1.5e-4)
# params.add('amp2', value = -3e-5, max = 0)
# params.add('c2', value = ca_reg_tweak[1])
# params.add('std2', min=0, value = 1.5e-4)
# params.add('amp3', value = -3e-5, max = 0)
# params.add('c3', value = ca_reg_tweak[2])
# params.add('std3', min=0, value = 1.5e-4)
# # Add parameters for the Voigt profile
# # params.add('m4_amplitude', value=-5e-5, max=0)
# # params.add('m4_center', value=ca_reg_tweak[3])
# # params.add('m4_sigma', value=1.5e-4)
# # params.add('m4_gamma', value=1.5e-4)
# # params.add('m4_skew', value =1.5e-4)
# # tweak lab wavelengths
# params['c1'].expr = f'{ca_reg_tweak[0]}'
# params['c2'].expr = f'{ca_reg_tweak[1]}'
# params['c3'].expr = f'{ca_reg_tweak[2]}'
# # params['m4_center'].expr = f'{ca_reg_tweak[3]}'
# params.add('beta', value=1.) # beta = (1+v/c)
# # Link the Lorentzian center to beta
# # params['m4_center'].expr = f'{ca_reg_tweak[3]} * beta'
# result = model.fit(norm_flux_stack[i], params,
# x = wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i])
# # weights = 1/flux_unc_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i])
# result_stack.append(result) # lmfit ModelResult objects
# dely = result.eval_uncertainty(sigma=3) # lmfit uncertainty
# dely_stack.append(dely)
# # Append the best_model, error and parameters stacks
# # could probably make a loop to the tune of for value in result etc...
# params_stack.extend([result.params['amp1'].value, result.params['c1'].value, result.params['std1'].value,
# result.params['amp2'].value, result.params['c2'].value, result.params['std2'].value,
# result.params['amp3'].value, result.params['c3'].value, result.params['std3'].value])
# # result.params['m4_amplitude'].value, result.params['m4_center'].value, result.params['m4_sigma'].value])
# beta_stack.append(result.params['beta'].value)
# # params_error_stack.extend([result.params['amp1'].stderr, result.params['c1'].stderr, result.params['std1'].stderr,
# # result.params['amp2'].stderr, result.params['c2'].stderr, result.params['std2'].stderr,
# # result.params['amp3'].stderr, result.params['c3'].stderr, result.params['std3'].stderr,
# # result.params['amp4'].stderr, result.params['c4'].stderr, result.params['std4'].stderr])
# best_model_stack.append(result.best_fit)
# mean_sigma_stack.append(np.mean(a =(result.params['std1'].value,
# result.params['std2'].value,
# result.params['std3'].value)))
# # result.params['m4_sigma'].value)))
# # print(result.fit_report())
# # print(f"Widths: {result.params['std1'].value:.8f}, {result.params['std2'].value:.8f}, {result.params['std3'].value:.8f}, {result.params['std4'].value:.8f}")
# # print(f"Mean Width: {np.mean(a= (result.params['std1'].value, result.params['std2'].value, result.params['std3'].value, result.params['std4'].value)):.8f}")
# params_arr = np.array(params_stack)
# # 0: amplitude
# # 1: centers
# # 2: sigmas
# amps = params_arr[0::3] # grab every amplitude
# amps_err = params_error_stack[0::3]
# centers = params_arr[1::3] # grab every center wavelength
# center_err = params_error_stack[1::3]
# sigmas = params_arr[2::3] # grab every sigma
# sigmas_err = params_error_stack[2::3]
# # print(centers)
# # Gaussian model fits for each source
# ca1_2_2614_fits = []
# ca1_2_2631_fits = []
# ca1_2_2657_fits = []
# fe1_2_2626_fits = []
# num_gauss = 3
# # # Getting the individual Gaussians into their own lists
# for i in range(len(standard_table)):
# ca1_2_2614_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[0::num_gauss][i], centers[0::num_gauss][i]*beta_stack[i], sigmas[0::num_gauss][i]))
# fe1_2_2626_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[1::num_gauss][i], centers[1::num_gauss][i]*beta_stack[i], sigmas[1::num_gauss][i]))
# ca1_2_2631_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[2::num_gauss][i], centers[2::num_gauss][i]*beta_stack[i], sigmas[2::num_gauss][i]))
# # ca1_2_2657_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[3::num_gauss][i], centers[3::num_gauss][i]*beta_stack[i], sigmas[3::num_gauss][i]))
# # Plotting each component Gaussian and the best model fit
# for i in range(len(standard_table)):
# # ca_reg_tweak = ca_reg_lines + np.mean(delt_diff_stack[i])
# # Create subplots with adjusted size ratios
# fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)
# ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], best_model_stack[i], c='red', label=r'Best Fit $\pm 3 \sigma$')
# ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], norm_flux_stack[i], c='black', alpha=0.5)
# # ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], result_stack[i].init_fit, '--', c='blue' , alpha=0.2, label='initial fit')
# ax1.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
# best_model_stack[i]-dely_stack[i], best_model_stack[i]+dely_stack[i],
# alpha=0.2, color = 'red' )
# # plot each componenet
# ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2614_fits[i], ls='-.',
# label=r'$\bf{Ca~I~2.2614~\mu m}$' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[0::num_gauss][i]*beta_stack[i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")
# # rf"$\pm {amps_err[0::num_gauss][i]:.3e}$"
# ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], fe1_2_2626_fits[i], ls='-.',
# label=r'$\bf{Fe~I~2.2626~\mu m}$' "\n" rf"A={amps[1::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[1::num_gauss][i]*beta_stack[i]:.5f}" "\n" rf"$\sigma$= {sigmas[1::num_gauss][i]:.3e}")
# ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2631_fits[i], ls='-.',
# label=r'$\bf{Ca~I~2.2631~\mu m}$' "\n" rf"A={amps[2::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[2::num_gauss][i]*beta_stack[i]:.5f}" "\n" rf"$\sigma$= {sigmas[2::num_gauss][i]:.3e}")
# # ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2657_fits[i], ls='-.', color='goldenrod',
# # label=r'$\bf{Ca~I~2.2657~\mu m}$' "\n" rf"A={amps[3::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[3::num_gauss][i]*beta_stack[i]:.5f}" "\n" rf"$\sigma$= {sigmas[3::num_gauss][i]:.3e}")
# # ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], best_model_stack[i]-ca1_2_2631_fits[i]-fe1_2_2626_fits[i]-ca1_2_2614_fits[i]+3, ls='-.', color='goldenrod',
# # label=r'$\bf{Ca~I~2.2657~\mu m}$' "\n" rf"A={result_stack[i].params['m4_amplitude'].value:.3e}" "\n" rf"$\mu$={result_stack[i].params['m4_center'].value:.5f}" "\n" rf"$\sigma$= {result_stack[i].params['m4_sigma'].value:.3e}")
# for j in range(len(ca_reg_lines)):
# ax1.axvline(ca_reg_lines[j]*beta_stack[i], ymin=0.05, ymax=0.95, color='black',alpha=0.4, ls='-.')
# ax1.axvline(ca_reg_tweak[j]*beta_stack[i], ymin=0.05, ymax=1., color='black', ls='-.')
# plt.suptitle(rf"Ca Region {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
# ax1.set_ylabel('Flux + Constant')
# ax1.legend(bbox_to_anchor=(1,1), frameon=False)
# ax1.set_title(r'Free $\sigma$ and $\beta$ offset, w/ different model on last line')
# # Plot residuals
# # ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
# # result_stack[i].plot_residuals(ax=ax2, datafmt='r.',markersize=3)
# ax2.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], result_stack[i].residual, 'r.', markersize=3)
# ax2.axhline(y=0, color='grey', alpha=0.5, linestyle='--')
# # ax2.set_ylim(-0.1,0.1)
# ax2.set_xlabel(r'Wavelength [$\mu$m]')
# ax2.set_ylabel('res')
# ax2.set_title(rf"$R^2$={result_stack[i].rsquared:.3f}, " r'$\beta \equiv 1+\frac{v}{c}= $' f"{beta_stack[i]:.8f}", fontsize=12)
# plt.show()
All Sources¶
# Size of 1 spectral resolution element
# IGRINS Spectral Resolution
spec_res = 0.00001
# Reduced and order-merged data filepath
# Laptop Path
data_path = "C:\\Users\\Savio\\Documents\\GitHub\\IGRINS-Spectra\\IGRINS_Merged"
# File path for figures to live in
fig_path = "C:\\Users\\Savio\\Documents\\GitHub\\IGRINS-SpectraIGRINS_figs\\standards_spectra"
# Create the folder if it doesn't exist
if not os.path.exists(fig_path):
os.makedirs(fig_path)
# Nicole's merged K-band spectra of some Taurus Standards
# merged_standard_files = glob.glob(data_path + "/merged_standards/m*.fits")
standard_table = pd.read_csv('C:\\Users\\Savio\\Documents\\IGRINS-Spectra\\standard_table_v3.txt', index_col=0) # csv of standards with file and Spectral Type, c/v TBA
proto_table = pd.read_csv('C:\\Users\\Savio\\Documents\\IGRINS-Spectra\\protostar_table.txt', index_col=0)
standards_path = standard_table['File']
standard_list = standard_table['File'].values
proto_path = proto_table['File']
proto_list = proto_table['File'].values
standard_stack = [] # initialize list of IGRINSSpectrum objects
proto_stack = []
num_orders = 25
order_len = 2048
for i in tqdm(range(len(standard_table))):
spec = IGRINSSpectrumList.read(standards_path[i])\
.trim_overlap(pivot=0.7)\
.remove_nans().stitch()
standard_stack.append(spec)
for i in tqdm(range(len(proto_table))):
spec = IGRINSSpectrumList.read(proto_path[i])\
.trim_overlap(pivot=0.7)\
.remove_nans().stitch()
proto_stack.append(spec)
# Determine the maximum length of flux arrays for the standards
max_flux_length = max(len(spec.flux.value) for spec in standard_stack)
max_wavelen_length = max(len(spec.wavelength.value) for spec in standard_stack)
max_snr_length = max(len(spec.snr.value) for spec in standard_stack)
max_unc_length = max(len(spec.uncertainty.array) for spec in standard_stack)
# Initialize stacks with NaN values
wavelen_stack = np.full((max_wavelen_length, len(standard_table)), np.nan)
flux_stack = np.full((max_flux_length, len(standard_table)), np.nan)
snr_stack = np.full((max_snr_length, len(standard_table)), np.nan)
flux_unc_stack = np.full((max_unc_length, len(standard_table)), np.nan)
for i in tqdm(range(len(standard_table))):
# convert angstrom to micron
wavelen = standard_stack[i].wavelength.value*1e-4
flux = standard_stack[i].flux.value
snr = standard_stack[i].snr.value
raw_flux_err = standard_stack[i].uncertainty.array
wavelen_mask = (wavelen > 2.1) & (wavelen < 2.39)
flux_mask = (flux > 0)
# flux_err_mask = (raw_flux_err < 0.1)
mask = (wavelen_mask) & (flux_mask) #& (flux_err_mask)
wavelen = wavelen[mask]
flux = flux[mask]
snr = snr[mask]
raw_flux_err = raw_flux_err[mask]
wavelen_stack[:len(wavelen),i] = wavelen
flux_stack[:len(flux),i] = flux
snr_stack[:len(snr),i] = snr
flux_unc_stack[:len(raw_flux_err),i] = raw_flux_err
# Directly query NIST to find line features in K-band
with warnings.catch_warnings(): # Ignore warnings
warnings.simplefilter('ignore')
lines_table = Nist.query(2.08*u.um,2.35*u.um,
linename = 'Na I, Sc I, Si I, Fe I, Al I, Mg I, Ca I, H I, Ti I',
energy_level_unit='eV',output_order='wavelength')
igrins_wav_cut = (lines_table['Observed'] > 2.08) & (lines_table['Observed'] < 2.35)
lines_table = lines_table[igrins_wav_cut]
# lines_table = pd.read_csv('lines_table.txt')
# Make masks for the table of all the lines just in case I want to peek at certain transitions/wavelengths
na1_mask = lines_table['Spectrum'] == 'Na I'
sc1_mask = lines_table['Spectrum'] == 'Sc I'
si1_mask = lines_table['Spectrum'] == 'Si I'
fe1_mask = lines_table['Spectrum'] == 'Fe I'
al1_mask = lines_table['Spectrum'] == 'Al I'
mg1_mask = lines_table['Spectrum'] == 'Mg I'
ca1_mask = lines_table['Spectrum'] == 'Ca I'
h1_mask = lines_table['Spectrum'] == 'H I'
ti1_mask = lines_table['Spectrum'] == 'Ti I'
# Just add all the masks to a list for the sake of my plotting a few cells down
mask_list = [na1_mask,sc1_mask,si1_mask,fe1_mask,al1_mask,mg1_mask,ca1_mask,h1_mask,ti1_mask]
color_list = ['purple', 'orange', 'green', 'blue', 'brown', 'crimson', 'olive', 'cyan', 'darkgreen']
na_reg_lines = sorted([lines_table[na1_mask][0]['Observed'],
lines_table[sc1_mask][19]['Observed'],
lines_table[si1_mask][2]['Observed'],
lines_table[sc1_mask][20]['Observed'],
lines_table[na1_mask][1]['Observed']])
na_reg_lines = np.array(na_reg_lines)
ca_reg_lines = list(lines_table[ca1_mask]['Observed'][:]) # list of Ca I lab wavelengths
ca_reg_lines.append(lines_table[fe1_mask]['Observed'][104])
ca_reg_lines = np.array(sorted(ca_reg_lines))
0%| | 0/23 [00:00<?, ?it/s]100%|██████████| 23/23 [00:07<00:00, 3.06it/s] 100%|██████████| 48/48 [00:15<00:00, 3.04it/s] 100%|██████████| 23/23 [00:00<00:00, 560.53it/s]
ca_reg_lines = list(lines_table[ca1_mask]['Observed'][:]) # list of Ca I lab wavelengths
ca_reg_lines.append(lines_table[fe1_mask]['Observed'][104])
ca_reg_lines = np.array(sorted(ca_reg_lines))
# Define the region for fitting
line_name = lines_table[ca1_mask][0]['Spectrum'] # Species
line_center = ca_reg_lines[1] # Wavelength
# from igrins_mod import local_continuum_fit
continuum_stack = []
continuum_err_stack = []
norm_flux_stack = []
norm_flux_err_stack = []
reg_idx_stack = []
# regions [left point,width]
regions = [(-1000,11), (-500,250), (190,51), (500,11) , (1000,11)]
# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
poly_deg = 3
for i in tqdm(range(len(standard_table))):
fig = plt.figure(figsize=(15,3))
# Perform the continuum fit
continuum, continuum_err, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
flux_stack[:,i],
flux_unc_stack[:,i],
poly_order = poly_deg,
line_center = line_center,
spec_res = spec_res,
regions = regions)
# Append indices to list of indices
reg_idx_stack.append(region_indices)
# Append to list of local continuum arrays
continuum_stack.append(continuum)
continuum_err_stack.append(continuum_err)
# Normalize flux by dividing the original flux by the continuum
norm_flux = flux_stack[region_indices[0][0]:region_indices[n][1], i] / continuum
norm_flux_stack.append(norm_flux)
# Plot the raw flux and continuum
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1], i], flux_stack[region_indices[0][0]:region_indices[n][1], i])
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1], i], continuum_stack[i])
# Highlight the regions and Ca I lines
for j in range(len(region_indices)):
plt.axvspan(wavelen_stack[region_indices[n - j][0], i], wavelen_stack[region_indices[n - j][1], i], color='red', alpha=0.2)
for k in range(len(ca_reg_lines)):
plt.axvline(ca_reg_lines[k], ymin=0.8, ymax=0.9, color='black')
plt.title(f"{line_name} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
plt.show()
0%| | 0/23 [00:00<?, ?it/s]
4%|▍ | 1/23 [00:00<00:06, 3.18it/s]
9%|▊ | 2/23 [00:00<00:05, 3.59it/s]
13%|█▎ | 3/23 [00:00<00:05, 3.61it/s]
17%|█▋ | 4/23 [00:01<00:05, 3.79it/s]
22%|██▏ | 5/23 [00:01<00:04, 3.84it/s]
26%|██▌ | 6/23 [00:01<00:04, 3.87it/s]
30%|███ | 7/23 [00:02<00:05, 2.67it/s]
35%|███▍ | 8/23 [00:02<00:05, 2.89it/s]
39%|███▉ | 9/23 [00:02<00:04, 3.17it/s]
43%|████▎ | 10/23 [00:03<00:03, 3.30it/s]
48%|████▊ | 11/23 [00:03<00:03, 3.53it/s]
52%|█████▏ | 12/23 [00:03<00:02, 3.67it/s]
57%|█████▋ | 13/23 [00:03<00:02, 3.86it/s]
61%|██████ | 14/23 [00:03<00:02, 3.90it/s]
65%|██████▌ | 15/23 [00:04<00:01, 4.07it/s]
70%|██████▉ | 16/23 [00:04<00:01, 4.04it/s]
74%|███████▍ | 17/23 [00:04<00:01, 3.09it/s]
78%|███████▊ | 18/23 [00:05<00:01, 3.23it/s]
83%|████████▎ | 19/23 [00:05<00:01, 3.46it/s]
87%|████████▋ | 20/23 [00:05<00:00, 3.64it/s]
91%|█████████▏| 21/23 [00:05<00:00, 3.80it/s]
96%|█████████▌| 22/23 [00:06<00:00, 3.83it/s]
100%|██████████| 23/23 [00:06<00:00, 3.57it/s] 100%|██████████| 23/23 [00:06<00:00, 3.57it/s]
Method 3¶
# Initialize storage lists
params_stack = []
beta_stack = []
params_error_stack = []
best_model_stack = []
result_stack = []
dely_stack = []
# For testing purposes
mean_sigma_stack = []
# ca_obs = ca_lines*(1+standard_table['v/c'][0]) # observed wavelengths (theoretically)
for i in range(len(standard_table)):
# Define initial parameters for Gaussian fitting
params = Parameters()
# Define Model
model = Model(ig.four_gaussian_beta, nan_policy='omit')
params.add('amp1', value = -1e-5, max = 0)
params.add('c1', value = ca_reg_tweak[0])
params.add('std1', min=0, value = 1.5e-4)
params.add('amp2', value = -1e-5, max = 0)
params.add('c2', value = ca_reg_tweak[1])
params.add('std2', min=0, value = 1.5e-4)
params.add('amp3', value = -1e-5, max = 0)
params.add('c3', value = ca_reg_tweak[2])
params.add('std3', min=0, value = 1.5e-4)
params.add('amp4', value = -1e-5, max = 0)
params.add('c4', value = ca_reg_tweak[3])
params.add('std4', value = 1.5e-4)
# tweak lab wavelengths
params['c1'].expr = f'{ca_reg_tweak[0]}'
params['c2'].expr = f'{ca_reg_tweak[1]}'
params['c3'].expr = f'{ca_reg_tweak[2]}'
params['c4'].expr = f'{ca_reg_tweak[3]}'
params.add('beta', value=1.) # beta = (1+v/c)
params['std2'].expr = 'std1' # constrain widths to be the same
params['std3'].expr = 'std1'
params['std4'].expr = 'std1'
result = model.fit(norm_flux_stack[i], params,
x = wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i])
# weights = 1/flux_unc_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i])
result_stack.append(result) # lmfit ModelResult objects
dely = result.eval_uncertainty(sigma=3) # lmfit uncertainty
dely_stack.append(dely)
# Append the best_model, error and parameters stacks
# could probably make a loop to the tune of for value in result etc...
params_stack.extend([result.params['amp1'].value, result.params['c1'].value, result.params['std1'].value,
result.params['amp2'].value, result.params['c2'].value, result.params['std2'].value,
result.params['amp3'].value, result.params['c3'].value, result.params['std3'].value,
result.params['amp4'].value, result.params['c4'].value, result.params['std4'].value])
beta_stack.append(result.params['beta'].value)
params_error_stack.extend([result.params['amp1'].stderr, result.params['c1'].stderr, result.params['std1'].stderr,
result.params['amp2'].stderr, result.params['c2'].stderr, result.params['std2'].stderr,
result.params['amp3'].stderr, result.params['c3'].stderr, result.params['std3'].stderr,
result.params['amp4'].stderr, result.params['c4'].stderr, result.params['std4'].stderr])
best_model_stack.append(result.best_fit)
mean_sigma_stack.append(np.mean(a=(result.params['std1'].value, result.params['std2'].value, result.params['std3'].value, result.params['std4'].value)))
# print(result.fit_report())
# print(f"Widths: {result.params['std1'].value:.8f}, {result.params['std2'].value:.8f}, {result.params['std3'].value:.8f}, {result.params['std4'].value:.8f}")
# print(f"Mean Width: {np.mean(a= (result.params['std1'].value, result.params['std2'].value, result.params['std3'].value, result.params['std4'].value)):.8f}")
params_arr = np.array(params_stack)
# 0: amplitude
# 1: centers
# 2: sigmas
amps = params_arr[0::3] # grab every amplitude
amps_err = params_error_stack[0::3]
centers = params_arr[1::3] # grab every center wavelength
center_err = params_error_stack[1::3]
sigmas = params_arr[2::3] # grab every sigma
sigmas_err = params_error_stack[2::3]
# print(centers)
# Gaussian model fits for each source
ca1_2_2614_fits = []
ca1_2_2631_fits = []
ca1_2_2657_fits = []
fe1_2_2626_fits = []
num_gauss = 4
# # Getting the individual Gaussians into their own lists
for i in range(len(standard_table)):
ca1_2_2614_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[0::num_gauss][i], centers[0::num_gauss][i]*beta_stack[i], sigmas[0::num_gauss][i]))
fe1_2_2626_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[1::num_gauss][i], centers[1::num_gauss][i]*beta_stack[i], sigmas[1::num_gauss][i]))
ca1_2_2631_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[2::num_gauss][i], centers[2::num_gauss][i]*beta_stack[i], sigmas[2::num_gauss][i]))
ca1_2_2657_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[3::num_gauss][i], centers[3::num_gauss][i]*beta_stack[i], sigmas[3::num_gauss][i]))
# Plotting each component Gaussian and the best model fit
for i in range(len(standard_table)):
# Create subplots with adjusted size ratios
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], best_model_stack[i], c='red', label=r'Best Fit $\pm 3 \sigma$')
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], norm_flux_stack[i], c='black', alpha=0.5)
# ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], result_stack[i].init_fit, '--', c='blue' , alpha=0.2, label='initial fit')
ax1.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
best_model_stack[i]-dely_stack[i], best_model_stack[i]+dely_stack[i],
alpha=0.2, color = 'red' )
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2614_fits[i], ls='-.',
label=r'$\bf{Ca~I~2.2614~\mu m}$' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[0::num_gauss][i]*beta_stack[i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")
# rf"$\pm {amps_err[0::num_gauss][i]:.3e}$"
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], fe1_2_2626_fits[i], ls='-.',
label=r'$\bf{Fe~I~2.2626~\mu m}$' "\n" rf"A={amps[1::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[1::num_gauss][i]*beta_stack[i]:.5f}" "\n" rf"$\sigma$= {sigmas[1::num_gauss][i]:.3e}")
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2631_fits[i], ls='-.',
label=r'$\bf{Ca~I~2.2631~\mu m}$' "\n" rf"A={amps[2::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[2::num_gauss][i]*beta_stack[i]:.5f}" "\n" rf"$\sigma$= {sigmas[2::num_gauss][i]:.3e}")
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2657_fits[i], ls='-.', color='goldenrod',
label=r'$\bf{Ca~I~2.2657~\mu m}$' "\n" rf"A={amps[3::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[3::num_gauss][i]*beta_stack[i]:.5f}" "\n" rf"$\sigma$= {sigmas[3::num_gauss][i]:.3e}")
for j in range(len(ca_reg_lines)):
ax1.axvline(ca_reg_lines[j], ymin=0.05, ymax=0.95,ls='--', color='black',alpha=0.5)
ax1.axvline(ca_reg_tweak[j]*beta_stack[i], ymin=0.05, ymax=0.95,ls='--', color='black')
plt.suptitle(rf"Ca Region {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
ax1.set_ylabel('Flux + Constant')
ax1.legend(bbox_to_anchor=(1,1))
ax1.set_title(r'Free $\sigma$ and yes $\beta$ offset')
# Calculate residuals
residuals = (norm_flux_stack[i] - best_model_stack[i])
# sqsum_res = np.sum(residuals**2)
# sqsum = np.sum((norm_flux_stack[i]-np.mean(norm_flux_stack[i])**2))
# R2 = 1-(sqsum_res/sqsum)
# Plot residuals
# ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
result_stack[i].plot_residuals(ax=ax2, datafmt='r.')
# ax2.axhline(y=0, color='grey', linestyle='--')
# ax2.set_ylim(-0.1,0.1)
ax2.set_xlabel(r'Wavelength [$\mu$m]')
ax2.set_ylabel('res')
ax2.set_title(rf"$R^2$={result_stack[i].rsquared:.3f}, " r'$\beta \equiv 1+\frac{v}{c}= $' f"{beta_stack[i]:.8f}", fontsize=12)
# Define the folder path
folder_path = 'Standards/ca_region'
# Create the folder if it doesn't exist
os.makedirs(folder_path, exist_ok=True)
# Save the plot to the file within the folder
plt.savefig(f'{folder_path}/{standard_table["Name"][i]}_linefit.png')
plt.show()
ew1_stack = []
ew2_stack = []
ew3_stack = []
ew4_stack = []
for i in range(len(standard_table)):
ew1 = np.trapz(np.abs(ca1_2_2614_fits[i]), wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew2 = np.trapz(np.abs(fe1_2_2626_fits[i]), wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew3 = np.trapz(np.abs(ca1_2_2631_fits[i]), wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
# ew4 = np.trapz(np.abs(1-ca1_2_2657_fits[i]), wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew1_stack.append(ew1)
ew2_stack.append(ew2)
ew3_stack.append(ew3)
# ew4_stack.append(ew4)
fig = plt.figure(figsize=(15,5))
# for i in range(len(standard_table)):
plt.plot(standard_table['Spectral_Type'], ew1_stack, ls='--', marker='o', c ='k', label = 'Ca I 2.2614')
plt.plot(standard_table['Spectral_Type'], ew2_stack, ls='--', marker='o', c ='red', label = 'Fe I 2.2626')
plt.plot(standard_table['Spectral_Type'], ew3_stack, ls='--', marker='o', c ='blue', label = 'Ca I 2.2631')
# plt.plot(standard_table['Spectral_Type'], ew4_stack, ls='--', marker='o', c ='green', label = 'Ca I 2.2657')
plt.xlabel('Source Spectral Type')
plt.ylabel(r'Equivalent Width [$\mu m$]')
# Define the folder path
folder_path = 'Standards/ca_region'
# Create the folder if it doesn't exist
os.makedirs(folder_path, exist_ok=True)
# Save the plot to the file within the folder
plt.savefig(f'{folder_path}/ca_region_{standard_table["Name"][i]}_eqwidth.png')
plt.legend(bbox_to_anchor=(1,1))
plt.show()
fig = plt.figure(figsize=(15,5))
# for i in range(len(standard_table)):
# plt.plot(standard_table['Spectral_Type'], np.array(ew1_stack)/np.array(ew1_stack), ls='--', marker='o', c ='k', label = 'Ca I 2.2614')
plt.plot(
standard_table['Spectral_Type'],
np.array(ew1_stack)/np.array(ew2_stack),
ls='--',
marker='o',
c ='red',
label = 'Ca I 2.2614/Fe I 2.2626'
)
plt.plot(
standard_table['Spectral_Type'],
np.array(ew1_stack)/np.array(ew3_stack),
ls='--',
marker='o',
c ='blue',
label = 'Ca I 2.2614/Ca I 2.2631'
)
# plt.plot(standard_table['Spectral_Type'], ew4_stack, ls='--', marker='o', c ='green', label = 'Ca I 2.2657')
plt.xlabel('Source Spectral Type')
plt.ylabel(r'Equivalent Width Ratios')
# Define the folder path
folder_path = 'Standards/ca_region'
# Create the folder if it doesn't exist
os.makedirs(folder_path, exist_ok=True)
# Save the plot to the file within the folder
plt.savefig(f'{folder_path}/ca_region_{standard_table["Name"][i]}_ew_ratios.png')
plt.legend(bbox_to_anchor=(1,1))
plt.show()